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INTRODUCTION 


I . 

Science  Applications,  Inc.  (SAI)  has  been  examining 
and  predicting  crater  ejecta  environments  for  both  high 
explosive  (HE)  and  nuclear  detonations  for  a number  of 
years.  Supported  by  the  Defense  Nuclear  Agency  (DNA), 

SAI  has  developed  a basic  ejecta  model,  recently  upgraded 
to  include  the  latest  theoretical  and  experimental  results 
(see  Reference  1).  The  work  reported  here  is  in  direct 
support  of  this  model  development  effort,  and  completes 
our  efforts  at  characterizing  the  ejecta  in-flight  size 
distributions  for  the  MIDDLE  GUST  series  of  high  explosive 
tests . 

The  SAI-developed  model  considers  the  entire  history 
of  crater  formation  and  ejecta  deposition.  Crater  forma- 
tion and  ejecta  throwout  are  modeled  and  trajectories  of 
representative  ejecta  are  calculated  to  provide  a descrip- 
tion of  the  airborne  ejecta  cloud.  The  cloud  is  followed  to 
impact,  and  a final  ejecta  blanket  thickness  is  obtained 
by  integrating  all  of  the  ejecta  impacting  the  surface. 

The  final  ejecta  blanket  is  generally  characterized  as 
having  two  distinctly  different  zones.  An  inner  zone 
which  can  be  described  as  having  a continuous  distribution 
of  ejecta  material  which  completely  covers  the  surface, 
and  an  outer  zone  where  impact  becomes  less  numerous,  the 
covering  is  discontinuous,  and  impact  assessments  are 
best  handled  in  a probabilistic  manner. 

The  work  reported  here  obtained  in-flight  ejecta 
size  distributions  from  photographic  records  of  the  MIDDLE 
GUST  events.  Combustion  products  (smoke)  and  large  amounts 
of  lofted  dust  tend  to  obscure  that  part  of  the  ejecta 
cloud  closest  to  detonation.  For  this  reason  these  results 
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are  more  properly  correlated  with  the  ejecta  found  in  the 
outer  impact  zone,  and  then  only  for  the  larger  ejecta 
as  the  limiting  resolution  of  the  photographic  recording 
system  will  not  permit  the  observation  of  ejecta  below  a 
certain  size  limit. 

The  following  two  sections  of  this  report  will 
outline  the  size  distribution  description  used  in  the 
ejecta  model,  will  review  the  supporting  data  for  it,  and 
will  then  report  on  the  size  distributions  measured  for 
the  MIDDLE  GUST  events.  The  report  concludes  with  a brief 
review  of  ejecta  ballistics  and  the  techniques  used  in 
deriving  dynamic  ejecta  parameters  such  as  velocity, 
position,  drug  coefficient  and  ejection  angle  from  simple 
photographic  records'. 


II.  EJECTA  SIZE  MODEL  AND  PAST  MEASUREMENTS 


The  SAI  ejecta  model  has  two  main  elements  --  the 
crater  source  and  a trajectory  calculation.  The  crater 
source  is  defined  by  a theoretical  velocity  distribution 
for  the  continuous  hydrodynamic  material  leaving  the 
crater,  by  empirical  models  for  ejecta  mass,  ejecta  size 
distribution  and  maximum  ejecta  size,  and  by  experimental 
observations  of  ejecta  velocity  (for  high  explosive  shots). 
The  trajectory  calculation  follows  the  ejecta  produced 
by  the  source  and,  using  a compressible  flow  aerodynamic 
drag  model,  tracks  them  through  a flowfield  defined  by 
the  VORDUM  (vortex  djjst  mode' 1 ) program.  Using  this  approach 
we  are  able  to  characterize  both  the  airborne  ejecta  cloud 
and  the  ground  deposition. 

The  work  reported  here  lias  been  directed  specifically 
at  that  element  of  the  model  which  covers  a description  of 
ejecta  size  distributions.  This  is  an  important  part  of  the 
ejecta  model,  and  is  currently  based  on  an  early  review 
(Reference  2)  of  many  different  sources  of  ejecta  --  con- 
ventional and  nuclear  detonations,  hypervelocity  impact 
studies,  meteor  impacts,  etc.  For  cohesive  media  such  as 
clay,  basalt  and  granite,  the  differential  mass  density 
distribution  is  assumed  to  be  represented  by 
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where  1 M/1  a is  the  ratio  of  the  di  fferential  mass 
associated  with  the  differential  ejecta  diameter 
am  is  the  diameter  of  the  largest  ejecta,  and  M 
total  ejecta  mass  in  megatons.  This  relationship 
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over  a range  of  ejecta  sizes  from  10  to  1 meter,  and 
Gault,  et  al  (Reference  3)  compiled  data  for  explosive 
and  impact  cratering  events  to  derive  a correlation  of 
mass  of  largest  ejecta  (M^)  with  as 

M = 6.6  x 10°  M 0-8 
m e 


(M  in  kg)  over  14  orders  of  magnitude  variation  in  M . 
Assuming  that  the  ejecta  can  be  represented  as  spherical 
(a  good  assumption  as  was  shown  in  Reference  4),  then 
M[n  = -npa^/6  and  the  above  relationship  would  imply  that 


)N  -3.5 
» — a a 
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where  3N/9a  is  the  ejecta  differential  number  density  as 
a function  of  size.  With  the  upper  limit  of  a , Equation  1 
is  assumed  valid  for  all  cratering  media  'xcept  those  specifi 
cally  identified  as  alluvial  or  sand  with  water  table  levels 
below  the  apparent  crater  depth.  For  the  latter  media, 
Reference  2 suggests  a more  accurate  relationship  of  the 


The  maximum  ejecta  diameter  for  such  soils  is  likely  to  be 
best  equated  with  maximum  in-situ  diameter  (rather  than  M(, ) , 
and  this  diameter  is  employed  as  the  upper  limit  in  Equation  4. 

The  experimental  data  upon  which  the  above  two  power 
laws  are  based  generally  falls  into  two  classes.  First,  there 
are  post-shot  observations  made  by  direct  retrieval  of 
impacting  ejecta  or  by  making  secondary  estimates  based 
upon  size  of  ejecta  impact  craters  or  upon  accumulation 
of  all  ejecta  debris.  Figure  1 is  a typical  summary  of 
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Figure  1 Size  Distributions  for  Five  NTS  Shots 


stub  data.  Second,  observations  have  been  made  by  in-si  Lu 
sampling  of  the  smaller  eject  a (dust)  which  are  lotted  and 
held  up  by  the  convective  forces  surrounding  the  rising 
fireball  (see  Figure  2).  Neither  of  these  techniques 
are  direct  observations  of  in-flight  characteristics  for 
larger  e.jecta,  and  a different  approach  was  sought  which 
would  provide  such  information  and  hopefully  validate  the 
use  of  the  above1  power  laws. 


Taken  from  Reference 


III. 


EJECTA  SIZING  FOR  THE  MIDDLE  GUST  EVENTS 


The  Denver  Research  Institute  (DRI),  under  contract 
to  DNA . has  been  tasked  with  obtaining  documentary  and 
technical  photography  of  almost  all  of  the  high  explosive 
tests  which  DNA  has  supported.  Of  specific  interest  are 
the  higher  framing  rate  records  which  have  been  utilized 
by  DRI  in  obtaining  fireball  and  dust  cloud  dimensional 
measurements  and  subsequent  drift  rates,  and  in  those 
records  which  were  obtained  with  exposure  times  and  spatial 
resolutions  sufficient  to  obtain  dynamic  parameters  for 
the  crater  ejecta.  The  latter  records  were  exposed  to 
permit  dimensional  measurements  with  enough  accuracy 
to  allow  the  estimation  of  ejecta  velocity,  ejection 
and  impact  angles,  and  of  drag  coefficient.  Coinci- 
dentally, these  records  provide  fairly  wide  field  "snap- 
shots" of  the1  airborne  ejecta  cloud  which  can  be  used 
to  quantize  in-flight  size  distributions  (see  Reference 
6 for  a detailed  description  of  the  photography,  and 
of  the  analysis  techniques  employed  by  DRI). 

Obtaining  statistically  meaningful  data  from  these 
records  is  simple  enough  in  concept,  but  difficult  in 
practice.  A direct  frontal  assault  would  imply  the 
utilization  of  a magnifying  lens  on  the  original  imagery 
; nd/or  controlled  enlargements  of  selected  frames  to 
make  manual  observations  of  the'  ejecta  sizes.  This  is 
obviously  a tedious  and  error-prone  way  to  proceed,  and 
it  was  determined  that  an  image  analysis  system  (specifically 


the  IMANCO  Quantimet)  could  perform  the  necessary  sizing 
and  analysis  task  (see  Reference  7)  with  a minimum  of 
observer  discomfort  and  with  increased  accuracy  and  through- 
put. An  initial  analysis  started  with  two  of  the  larger 
shots  (20  and  100  ton)  of  the  MIDDLE  GUST  high  explosive 
test  series  (see  References  4 and  8)  and  the  work  reported 
here  continues  with  ejecta  sizing  for  other,  lower  yield, 
events  of  that  same  series. 

The  events  analyzed  in  References  4 and  8 were 
MIDDLE  GUST  I , a 20  ton  half-buried  sphere  of  TNT,  and 
MIDDLE  GUST  III,  a 100  ton  tangent-above  sphere  of  TNT. 
Qualitatively,  these  two  shots  yielded  in-flight  sizing 
data,  shown  in  Figure  2,  which  was  similar.  The  a ' ° 
power  law  curve  superimposed  on  these  measured  curves 
might,  at  first  glance,  imply  that  the  use  of  such  a 
power  law  in  the  ejecta  model  is  questionable.  To 
resolve  this,  the  ejecta  model  was  exercized  lor  the  100 
ton  event  and  the  airborne  ejecta  field  at  5 seconds  was 
mapped  in  detail.  That  portion  of  the  ejecta  field  which 
coincided  with  the  f ield-of-view  of  DRI's  recording  cameras 
was  isolated  and  numerically  summed  to  create  a size 
distribution  curve  which  was  directly  comparable  to  the 
measured  data.  These  results  are  shown  in  Figure  1 
where  it  can  be  seen  that  there  is  reasonable  agreement 
between  the  prediction  and  the  experimental  observations. 

It  should  be  noted  that  at  least  some  of  the  flattening 
of  the  experimental  curve,  for  smaller  ejecta,  is  due  to 
resolution  problems  associated  with  the  recording  systems 
(film  and  camera)  and  with  the  image  analysis  system  (see 
Reference  4). 

With  such  general  validations  of  the  power  law 
size  distribution  in  the  ejecta  model,  it  was  thought 
appropriate  to  extend  such  experimental  observations  to 
other  geologies.  The  calibration  shots  for  the  MIDDLE 
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GUST  scries  were  1000  pound  TNT  detonations  which  were 
set  off  over  varying  sub-surface  geologies.  The  overall 
geology  was  that  for  events  I and  III  except  that  pre-shot 
excavating  permitted  detonations  over  changing  sub-surface 
media  as  is  shown  in  Figure  5.  The  four  shots  illustrated 
in  this  figure  were  selected  for  analysis  because  they  had 
moderately  good  photographic  records  and  because  the  only 
difference  between  the  shots  was  sub-surface  geology  , i . e . , 
all  four  shots  were  1000  pound  half-buried  spheres  of  center- 
initiated  TNT. 

As  for  shots  I and  III,  the  calibration  shots  were 
recorded  with  DRI's  70mm  Hulcher  cameras  with  slight 
changes  in  range- to-detonat ion  and  in  the  focal  length 
of  the  lenses  used.  The  pertinent  recording  parameters 
are  summarized  below,  but  all  appropriate  geometric 
corrections  have  been  included  in  the  results  presented. 


Table  1 


Pertinent  Recording  Parameters  for 
MIDDLK  GUST  Calibration  Events 

Event  Camera  Range  Lons  Focal  Charge  Position 


( ft) 

Length  (mm) 

from  Original 
Ground  Level 

3a 

■147 

1 50 

0 

G 

670 

1 50 

4 

8 

75  1 

150 

9 

9 

757 

150 

16 

For  all  four  shots,  the  airborne  ejecta  field  was  quantized 
at  2,  1 and  6 seconds  after  detonation.  On  the  3a  event 

the  recording  cameras  were  considerably  closer  to  the 
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detonation  than  for  the  other  shots  and  therefore  by  G 
seconds  after  detonation  much  of  the  e.jeeta  field  has 
exceeded  the  recording  cameras  f i e ld-of- v iew  and  no  data 
were  acquired.  Full  field  measurements  were  made,  either 
Fast  or  West  of  the  detonation,  and  several  symmetry  checks 
were  made  at  two  seconds  after  detonation.  For  all  <»1  these 
shots,  the  ejecta  field  had  not  exceeded  recording  bounds 
so  that  measurements  generally  included  all  ejecta,  on  one 
side,  up  to  the  edge  of  the  fireball/dust  cloud. 

Figures  G through  9 are  typical  of  the  data 
which  were  acquired  in  these  sizing  measurements.  For 
clarity,  the  data  points  have  been  eliminated,  and  only 
the  "eyeball"  fits  are  shown  to  indicate  general  trends. 

(Data  for  these  plots  are  tabulated  in  Appendix  A.)  A 
review  of  these  figures  will  show  that  shots  da,  G and  9 
have  very  similar  size  distribution  curves,  that  the  total 
ejecta  field  diminishes  in  time*,  and  that  the  larger  ejecta 
tend  to  drop  out  preferentially.  Also,  it.  would  appear 
that  shot  9,  over  a shale  geology,  tended  to  produce  signifi- 
cantly more  ejecta  than  either  ol  tin1  other  two  shots. 

Shot  8 appears  to  differ  from  the  others  in  only  one 
respect,  that  its  distribution  curves  appear  to  be  flatter 
(i.e.,  more  larger  ejecta,  or  fewer  small  ones).  The 
reason  for  this  is  not  apparent,  but  it  was  set  off  over 
a clay  geology  which  might  have  biased  the  distribution. 

Figures  10  and  11  take  the  two  and  four 
second  curves  from  the  preceding  four  figures,  and  combine 
them  for  direct  comparison.  Here  the  increased  ejecta 
count  for  shot  9 is  readily  seen,  particularly  at  two  seconds 
where  an  apparent  factor  of  two  in  total  count  is  readily 
observed.  The  differential  has  dropped  by  four  seconds, 
but  shot  9 is  still  on  the  upper  bound  of  the  envelope. 

Also  seen  in  these  figures  is  the  apparent  flattening  of 
the  distribution  lor  shot  8. 


Integral  Number  Observed 
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figure  10 


Ejecta  Area 


(cm2) 


Comparison  of  Ejecta  Sizing  at  2 Seconds  Alter  Detonation 
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Figures  12,  13,  and  14  compare  ejecta  field 

measurements  made  Fast  and  West  of  shots  6,  8 and  9,  two 

* 

seconds  alter  detonation.  It  can  be  seen  from  these 
figures  that  the  events  were  indeed  highly  symmetrical 
and  that  there  is  no  apparent  biasing  of  the  data  in  either 
of  the  viewing  directions. 

In  Figure  15  we  combine  size  distribution  curves 
for  the  four  calibration  shots  with  data  from  the  20  and 
100  ton  events.  The  results  for  the  calibration  shots  and 
the  20  ton  shot  should  be  directly  comparable  as  these  were 
all  half-buried,  center-initiated,  spherical  charges,  but 
tin1  100  ton  event  was  tangent-above  and  therefore  not 
directly  comparable. 


A similar  comparison  could  not  be  made  for  shot  3a  because 
a different  film,  which  was  underexposed,  had  been  used  on 
one  side.  It  should  also  be  noted  that  shots  G,  8 and  9 
were  recorded  on  Shellburst  film,  while  3a  was  recorded  on 
tri-x  panchromatic.  Data  presented  here  for  the  proceeding 
three  shots  are  therefore  directly  comparable,  while  that 
for  3a  might  be  slightly  biased  because  tri-x  does  not 
quite  have  the  spatial  resolution  capability  of  Shellburst. 


:i  . E.jec  I a Di  amet  < 


( cm ) 


i^urc  15  Comparison  of  L.jeota  Si/.inp:  for  MIDDLE 
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IV.  DETERMINATION  OF  EJECTA  DYNAMICS  FROM 
PHOTOGRAPHICALLY  RECORDED  TRAJECTORIES 


One  o 1 the  main  objectives  of  the  DRI  photographic 
projects  is  to  obtain  data  on  ejecta  dynamics.  Originally, 
these  investigations  were  limited  to  in-flight  and  impact 
parameters,  but  as  the  SAI  ejecta  model  was  developed  an 
additional  need  for  quantitative  ejection  parameters 
evolved.  To  obtain  the  desired  data  DRI  has  developed 
an  analysis  program  which  uses  trajectory  coordinates 
obtained  from  digitized  photographic  records  (Reference  9). 

Ejection  velocities  initially  reported  by  DRI  were 
inconsistent  with  values  being  used  in  the  ejecta  model 
and  implied  impact  ranges  which  were  considered  too  long. 

A critical  review  of  the  DRI  analysis  procedure  revealed 
several  defects,  and  their  approach  has  recently  been 
improved  to  provide  more  reasonable  and  accurate  ejection 
parameters  (Reference  6).  However,  the  ejecta  are  expelled 
from  the  crater  at  some  azimuthal  angle  ( i . e. , toward  or 
away  from  the  recording  camera)  and  determination  of  this 
angle  from  the  recorded  trajectory  remains  uncertain  at 
best.  SAI  was  asked  to  review  trajectory  analysis  tech- 
niques and  to  determine  an  optimum  approach  (i I any)  to 
measure  this  angle  and  therefore  to  improve  the  accuracy 
of  estimated  dynamic  parameters. 

The  coupled  differential  equations  of  motion  which 
describe  (in  the  trajectory  plane)  ballistic  trajectories 
lor  ejecta  are  given  by 

x = d'x/dt'  = “ y/(  dx  /d  t ) + ( dy  /d  t.  )2  dx/dt 


1 


where 


y°  = ci-y/dt.  ’ = v/cdx/dt)2+  (dy/dt)2  dy/dt 


a ml 


x.y  = horizontal  and  vertical  trajectory  coon 
g = gravity 

K = drag  force  = ( t'^A/nt 


P = 

CD  = 
A = 

m = 


a tmospher i c dens i ty 
drag  coefficient 
pro  .joe  t ed  ejecta  area 
eject' a mass. 


These  equations  may  be  used  to  analyze  ejecta 
trajectories  by  assuming  no  ablation  or  agglomeration 
(or  the  ejecta,  and  by  assuming  some  simple  model  for 
CD  (o.g.  constant).  Note  that  in  the  absence  ot  drag 
these  equations  reduce  to  those  describing  a simple 
parabolic  trajectory,  and  near  trajectory  peak  these 
equations  again  reduce  to  a simpler  form.  Near  peak, 
y is  maximum  and  dy/dt  0.  Then 


d ' x 
dt  2 


du 
d t 


K j d \ 
2 \dt 


K dx 
2 u dt 


Hew  r i t i ng  wo  have 


du  K , 

TT  = -5  dx 


a nd  i n t egra t i ng , 


u(\2  ) = u(  xt  ) exp  - r,  ( x,,  - Xj  ) 
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Near  the  trajectory  peak,  then,  the  horizontal  velocity 
component  is  degraded  exponentially  according  to  Equation  9. 
Note  that  the  above  equations  can  at  most  yield  K,  and 
that  a value  for  C'^  will  depend  on  the  accuracy  with  which 
ejecta  dimensions  may  be  measured. 

To  correctly  interpret  photographically  recorded 
ejecta  trajectories,  specifically  to  obtain  ejection,  impact 
and  dynamic  parameters,  requires  a detailed  understanding 
of  imaging  optics  and  the  geometric  relations  of  the  image 
plane,  the  plane-of-view , and  the  trajectory  plane.  Figures 
16  and  17  are  schematic  representations  of  the  image 
plane  (film)  and  the  plane-of-view  (note  that  we  have  here 
corrected  the  left-to-right  and  top- to-bot tom  inversion 
created  in  the  image  plane).  The  image  (film)  plane  is 
readily  understood,  but  we  now  define  the  plane-of-view 
as  that  plane  perpendicular  to  the  camera's  optic  axis 
and  passing  through  the  detonation  point  (GZ).  This  is 
the  simplest  plane  to  make  measurements  in,  and  the  distance 
from  camera  to  plane,  along  the  optic  axis,  is  called  the 
Range  on  Optic  Axis  ( ROA ) . 

In  general,  the  distance  from  camera  to  GZ  (the 
range,  R)  is  known,  as  is  t h€>  camera  aiming  in  terms  of 
elevation  and  azimuthal  angles  (El  and  Az)  with  respect 
to  GZ . R and  ROA  are  then  related  with  the  simple 
expressi on 

R = ROA  (1  + tan2  El  + tan2  Az)  . (10 

Also,  measurements  on  the  film  (d)  art?  directly  related 
to  distances  in  the  plane  of  vie-w  (D)  by 


Figure  16  Outline  of  Digitized  Ejecta  Trajectory 


where  m,  the'  magnification  factor,  is  defined  by 


m = ROA/f,  (12) 

and  f is  the  local  length  of  the  camera  lens.  Note  also 
that  A/,  and  El  may  be  readily  obtained  from  film  measurements 
of  the  x and  y displacements  of  GZ  from  the  optic  axis  as 


tan  El  = y(GZ)/f 

a n d 

tan  Az  = x(GZ)/f. 


(13) 

(ID 


ROA  is  then  calculated  with  Equation  lc  . 


It  is  now  noted  that  ejecta  measurements  in  the 
film  plane,  x(e)  and  y(e),  can  be  projected  onto  the 
plane-of-view  and  that  azimuthal  and  elevation  angles 
may  be  ascribed  to  the  ejecta  with  respect  to  the  optic 
axis  and  the  direction  vector  to  GZ . The  latter  is 
essentially  horizontal,  and  we  use  this  fact  and  the 
derived  azimuth  and  elevation  angles  to  determine  ejecta 
location  in  the  trajectory  plane. 


Referring  to  Figure  16,  the  center  "+"  is  meant 
to  represent  the  optic  axis  oi  the  recording  system.  The 
elevation  and  azimuth  angles  of  GZ  with  respect  to  the  optic 
a x i s a re  de t i ned  by 


and 


El(GZ)  = tan  1 

y(  GZ) 

t 

Az(GZ)  = tan  1 

x ( GZ  ) 
f 

(15) 


(16) 


Similarly,  the  El  and  Az  of  the  ejecta  with  respect  to 
the  optic  axis  are  defined  by 


F 1(e)  = 

tan"1 

(17) 

Az(e) 

tan’1 

( 1«) 

desired  angles 

0 and  P are  then  simply  given  by 

o = i; l ( v 

) - i;  L<  GZ ) 

(19) 

<P  = Az(  e 

■)  - Az(GZ). 

(20) 

This  is  illustrated  in  Figure  IS  where  the  azimuth 
and  elevation  angles  have  been  labeled  ip  and  0,  respectively, 
and  the  ejecta  plane  has  been  rotated  by  an  angle  3 with 
respect  to  a vertical  plane  perpendicular  to  the  line 
between  camera  and  GZ . The  angle  ft  i t se 1 I is  not  of 
interest,  but  its  knowledge  is  necessary  to  correctly 
calculate  the  x,y  position  of  ejecta  in  the  trajectory 
p 1 a ne . 

With  the  aid  of  Figure  18,  the  following  relationships 
are  readily  derived  to  determine  ejecta  position  in  the 
trajectory  plane  as  a function  of  R,  0,  ,,  and  ft: 


R s i n |_ 
>s(  * . ) 


(21  ) 


- R tan  J3 cos  ft  . , 

y ~ cos  (ft  + P)  * 

Cnrren  t Anal  ysis  T<  ‘chn  i gut  * 

Details  ol  the  analysis  scheme  currently  used  to 
derive  dynamic  parameters  from  digitized  imagery  art' 
found  in  Reference  (5,  This  approach  may  be  summarized 
in  tht'  following  four  steps: 


Figure  18  Outline  ot  Camera,  GZ  and  Trajectory  Plane 


1.  Assume  an  azimuthal  ejection  angle  (3)  for 
the  trajectory  to  be  analyzed. 

2.  Near  trajectory  peak,  use  Equation  9 to 
estimate  K by  fitting  the  trajectory  points 
with  a second  order  polynominal. 

3.  Use  this  3 and  K to  extrapolate  first  data 
points  toward  origin. 

4.  If  resultant  trajectory  originates  at  pre- 
determined location  (e.g.,  final  crater  radius) 
assume  correct  fit  in  and  K have  be  n found.  II 
not,  adjust  value  for  6 and  repeat  above  stops. 

This  approach  has  several  shortcomings  which  became  obvious 
during  an  early  review  of  the  problem.  The  derivation  of 
K is  biased  by  an  a-priori  second  order  polynominal  lit, 
and  t ho  value  derived  can  be  dependent  on  the  selected 
spacing  between  data  points.  Most  of  the  acquired  digital 
data  is  never  really  used  in  an  optimum  fashion  and 
forcing  the  trajectory  toward  a preconceived  origin  also 
seems  faulty.  It  was  felt  that  the  above  approach  could 
be  improved  upon,  and  that  perhaps  accurate  trajectory 
data  could  provide  valuable  insights  to  cratering  mechanisms 
and  therefore  the  accuracy  and  validity  of  existing 
cratering  models. 


An  Improved  Look  at  Kj< •da  Dynamics 

Realizing  the  deficiencies  of  the  above  approach, 

SAI  sought,  to  obtain  data  of  improved  accuracy  and  validity 
by  utilizing  a direct  approach  to  fit  al 1 acquired  trajectory 
data  in  a least  squares  test.  The  approach  can  be 
summarized  as  fol lows: 

Make  initial  estimates  of  the  dynamic  parameters 
sought  (including  3),  and  calculate  an  ejecta 


trajectory  through  the  coupled  equations 
ol  motion.  Compare  calculated  and  experi- 
mental trajectory  points  and  minimize,  in  a 
least  squares  sense,  the  difference  between 
the  two  sets  of  information  by  adjusting  the 
va 1 ues  of  the  dynamic  parameters  used  in 
calculating  the  trajectory. 

In  concept  the  above  approach  is  simple.  In  practice, 
because  of  the  coupled  differential  equations  of  motion 
(Equations  5 &,  (1 ) and  because  ol  the  trignometric  coordinate 
transformation  (Equations  21  x 22).  it  is  considerably  more 
difficult.  Development  of  a detailed  computer  analysis 
technique  to  perform  the  desired  calculations  would  have 
been  extremely  costly  and  time-consuming,  and  therefore 
beyond  the  scope  of  the  demonstration  sought.  This  rough 
scoping  operation  was  only  made  possible  by  utilizing  a 
proprietary  computer  language  called  PROSE,  which  is  a 
higher  order  language  capable  of  performing  direct 
calculus  calculations  (e.g.,  derivatives  are  explicitly 
calculated  and  carried  rather  than  the  digital  approximations 
which  are  so  familiar  to  most  FORTRAN  users).  Using  PROSE, 
we  were  able  to  examine  the  estimation  of  ejecta  dynamics 
from  the  digitized  trajectory  data,  and  to  make  some 
observations  on  the  validity  of  the  data  currently 
generated . 

As  mentioned  above,  the  criterion  for  best  estimate 
in  this  approach  is  to  minimize  the  difference  between 
measured  and  calculated  trajectory  points.  Letting  X(t.), 
Y(t.)  and  x ( t . ) , v(t)  bo  the  respective  measured  and 

i i l 

calculated  trajectory  points  at  times  t.,  we  seek  to 
minimize  the  sum  of  the  N straight  line  distances  between 
measured  and  calculated  trajectory  points: 


SL'MSQR 


Sum  of  squares  = 


N 

V 

i = l 


X(ti)-x(t.)] 2 + |Y(ti)-y(t.) 


(23) 


The  problem  was  formulated  in  PROSE  (see  Table  2 for 
typical  listing),  "ideal"  data  was  generated  through  the 
equations  of  motion  and  the  code  was  then  asked  to 
calculate  dynamic  parameters  from  selected  portions  of 
the  "ideal"  data. 

For  testing  purposes,  the  "ideal"  data  was  limited 
to  two  cases;  trajectories  initiated  at  a 15°  elevation 
angle  and  a velocity  of  200  feet  per  second,  both  with 
and  without  a drag  force.  The  case  with  a drag  force 
assumed  a "typical”  chunk  of  cohesive  dirt  with  a 30  cm 
diameter  and  0.6  so  that  K/2  = 9.56  ,\  10-ti(cgs).  The 

trajectories  generated  for  these  two  cases  are  shown  in 
Tables  3 and  1.  As  an  approximation  to  the  situation 
which  would  be  found  in  the  normal  analysis  procedure, 
these  two  trajectories  were  picked  up  at  three  seconds 
(i.e. , just  before  peak)  and  data  sets  of  from  10  to  20 
points  were  created.  The  PROSE- developed  code  (e.g., 

Table  2)  was  then  asked  to  calculate  some  or  all  of  the 
dynamic  parameters  at  3 seconds  which  could  be  used  as 
initial  conditions  to  match  the  "ideal"  data  through 
the  equations  of  motion.  Note  that  implicit  in  this 
approach  is  the  concept  that  given  the  initial  conditions 
at  3 seconds  which  can  be  used  to  duplicate  the  measured 
data,  then  these  same  conditions  can  be  used  to  back  the 
trajectory  toward  the  origin  to  determine  ejection 
parameters.  The  parameters  to  be  matched  in  these  two 
cases  are  listed  below  in  Table  5. 
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Table  2 


Typical  PROSE  Program  Listing  for 
Ejecta  Trajectory  Analysis 


: L 

131  ^ ED  GRAM  . 3 'JESS 

231  ALLOT  T^C 1 7) , XDC 1 7 ) , Y3t 1 7 ) ,30DC 6) , A3C 1 7 ) , EL ; 1 7 ) 

32  1 3 TAD  DATA 

42 1 SELECT  TTY 

52  1 30D  = . DATA  C 1 , 1 , 1 , 1 , 1 , 5 1 

63  1 IMPOSE  . SE  T 

73  1 FI  MD  \ Y 1 , - T A I M .TO  IP  WITH  30  OODS  • 1 

62  1 DOM?  OL0  3AL- 
2 3 1 E'l D 

133  1 MODEL  . T E I 3 
1131  T=T3 ; 1 ) 

1221  X=X3 
1321  Y=Y3 
1421  0=03 
152  1 V=V3 

1631  10 1 T!  ATE  JAOI3I3  “33  . D I E F TJOATIOOS 

17  3 1 ODOT/o,  o/X,  VDOT/ V , V/Y  IF  T S TE 0 OT 

1=21  S DM  S 0 2 = 1 

1331  FOE  1=2  TO  0 DO 

233  1 ’ ,XJ  T I L T 3E  T"5  ; I) 

212  1 ! 'IT  E 33  ATE  . D I rF 

223  1 5I0A T = . SI ODC A U I ) ) 

2 32  1 COS  3-=  . C05DC  3“T  A<- AT  < I ) > 

243  1 X3C  I > = 310 42/ CO  53° 

>533  Y?C I »=. TAODC ELC I ) ) »< 1+5 IOAT*. SI01  3ETA>/CC3 

2 '■'21  3 'JOS  3 2=  50X333  «■(  '<■=  C I)-X)  ♦ * 2 ♦ ; Y-3  C I ) -Y>  >*2 

2731  EE? FAT 

2521  EOD 

>23  1 MODEL  . DI17- 

332  1 "EL  = . lETt  U»D<-  V»V> 

313>  ODOT  = - ^ *"EL 
3231  ODOT=- r*0» VEL- GPAV 
3331  EOD 

343  1 CO  1T20LLE2  .SET  FOE  O-EA 

3 53  1 ADJ  OS  T=  2 
3551  DELTA  = 1 E-4 
362  1 C ) InE  >.GF  = 1 


Table  3 


Case  1 ( = 0.6)  Ideal  Trajectory  Used  for  Data  Simulation 


time 

(sec) 

dx/dt 
( cm/ sec ) 

dy/dt 
( cm/sec ) 

X 

(ft) 

y 

(ft) 

3 . 222E-2  1 

4.31  1 E+  3 3 

4.311 E + 3 3 

2 . 222E-2 1 

2.222E-21 

2.  523E-3  1 

4.  2 4 9 E + 3 3 

4.33  6E  + 3 3 

3.511 E+21 

3 . 4 1 5E+2 1 

5.  332E-2  1 

4. 1 92E+33 

3. 73  3E  + 3 3 

6.  9 74 E+  2 1 

6.  52 3 E+2  1 

7 . 52  <3  E - 2 1 

4. 1 37E+33 

3.  4!  6E  + 3 3 

1 . 3 3 9E  + 32 

9. 539E+2  1 

1 . 333E+23 

4.  3 3 5E+  3 3 

3 . 1 3 3 E * 3 3 

1 . 376E+22 

1 . 222E+22 

1 . 253E  + 33 

4.  336E+33 

2.  3 4 9 E+3  3 

1 . 739E+22 

1 . 46  5E  + 3 2 

1 . 533E+3  i 

3.  9 9 3 E +3  3 

2.  57  2 E +-3  3 

2 . 33  9E  + 22 

1 . 63 : E + 3 2 

1 . 7 5 3 E + 3 3 

3.  9 46E+3  3 

2.  33  3 E -t-3  3 

2 . 3 64 E+2 2 

1 .3  ■ j E + 2 2 

2. 3 3. 3 E+2  3 

3.  93  4E  + 3 3 

2 . 3 3 1 E + 3 3 

2.63  6E  + 22 

2. 266E+32 

2 . 2 53  E + 3 2 

3 . 3 6 3 E + 2 3 

1 . 767E  + 3 3 

3. 335E+32 

2 . 223E+22 

2. 533E+33 

3 . 5 2 5 E 3 i 

1 . 53  5E  +-3  3 

3.322E+32 

2 . 3 5 7 E + 2 2 

2.  753E+32 

3.  733  £+2 3 

1 . 247E+3 3 

3 . 632E+22 

2 . 4 72E+  2 2 

3. 332E+33 

i.  753E+33 

9.  93  3 E -t-3  2 

3 . 942E  + 32 

2. 563E+32- 

3. 252 E +3 3 

j . 7 1 3 - 3 3 

7 . 3 7 7 E + 3 2 

4.243E+22 

2 . 6 3 4 E+  3 2 

3 . 533E+33 

3.  6 j 5E+3  3 

4.  -3  72  E+3 2 

4. 552E+32 

~ • 654E+22 - 

3 . 753E+32 

3.  65  3 E+2 

2. 3 3 7 E+  3 2 

4 . 5 53  E+  2 2 

2.71 5E+22 

4. 332E+33 

3. 62 1 E+33 

-7. 496E+33 

5 . 1 5 1 E + 2 2 

2. 72+.E+32- 

4.  2 53E+33 

3.  593  E+3  3 

-2.51 5~+32 

5 . 4 4 7 E+  2 2 

2 . 7 1 4E+32 

A.  53  2 E +3  3 

3.  •.  9 £+3 

-4. 9357+32 

5.  742E  + 2 2 

2.  6 2 4 E + 2 2 — 

4. 753E+33 

■ . 29E+33 

- 7.  334E+3  2 

6 . 2 3 1 E + <.  j 

. . 63  4 E + 3 2 

5. 333E+33 

3.  493  E+3 

- 9 . 7 1 2 E + 3 2 

6.31 9E+2. 

2. 565E+22 

5.  2 53E  + 33 

3. 463E+3 3 

- 1 . 23  7E  +3  3 

6 . 62  5E  + 2 2 

2.476E+22 

5. 533E+33 

3.  437E  + 3 3 

- 1 • 443  E + 2 3 

2. 367E+32 - 

5. 753E+33 

3. 43  7 E +3  3 

-1.671 E + 3 3 

7. 1 69E+22 

2.243E+32 

6.333E+32 

3.  3 7 5 E+  3 3 

- 1 . 933E  + 33 

7.447E+2' 

2.2  9«iE  + 2 2 - 

6. 253E+23 

3. 344E+23 

- 2.  1 27  E + 3 3 

7. 723E+22 

1 . 9 2 9 E + 2 2 

6. 523E+33 

3.  3 1 2E*-2  3 

-2. 353E+33 

7 . 9 9 6E  + 3 . 

1 . 74 6 E+2  2 

6.  7 5 3 E + 3 3 

3.  232 E+2 3 

-2.  -571  E + 2 3 

J.266E+22 

1 . 5 4 5 E + 2 2 

7.  2 3 3 E +3  2 

3.  2 47 E* 23 

- 2.  73  9 E +2  3 

3. 534E+22 

1.525E+22  “ 

7.  2 53  E + 3 3 

3.  3 1 3 E + 3 3 

- 3 . 3 2 4 E + 3 3 

■ . 7 9 9 E + 2 : 

1 . J 2 i E + * 2 

7. 53 3 E + 3 3 

3.  1 7 9E  + 33 

-3. 216 E +33 

9 . 2 6 1 E + 2 _ 

j . 3 3 3 E + 2 1 

Table  4 


Case  2 (No  Drag)  Ideal  Trajectory  Used  for  Data  Simulation 


time 

dx/dt 

dy  / dt 

X 

y 

(sec ) 

( cm/ sec ) 

( cm/ sec ) 

(ft) 

(ft) 

3 . 3 3 3 E-  3 1 

4.311 

E + 3 3 

a . 

31  1 E + 3 3 

3 . 3 3 5 E - 3 1 

5.5-52-. 1 

2.  533E-3  1 

4.  3 1 1 

E+3  3 

4. 

365E  + 3 3 

3 . D 3 6E  + 3 1 

3 . 43  7E+5  1 

5. 333E-3 1 

4.311 

E + 3 3 

3. 

323E+33 

7. 3 71  E+3 1 

6. 6 75E+ . 1 

7.  5 33E-3  l 

4.311 

£♦■33 

3. 

57  5E+  3 3 

1 . 61  E+3  : 

9 . 7 3 3 E + . 1 

1 . 333E+33 

4.311 

E + 3 3 

3. 

3 33  E +■  3 3 

1 . 4 1 4 E + - : 

1 . 2 5 4 E + 5 . 

1 . 253E+33 

4.311 

E + 3 3 

3. 

J 3 5 E+  3 3 

1 . 7 6 3 E + 3 2 

1.51 7E+52 

1 .5 33E+33 

4.311 

E + 33 

2 . 

i43E  + 3 3 

2 . 1 2 1 E+J2 

1 . 7 6 1 E + . 2 

1 . 753E+33 

4.311 

E + 33 

2 

5 94E+33 

2.47aE+32 

1 • J ; 2 5 ♦c 

2. 333  E + 3 3 

4.  3 1 1 

E + 3 3 

2 

3 4 9 E + 3 3 

2 . 3 2 j E - 2 

2 . 1 ; 7 E + <_  _ 

2. 253E  + 33 

4.  3 1 1 

E + 3 3 

2 

1 3 4 E + 3 3 

3. 1 32E+32 

2 . 3 fc  + E + ..  - 

2. 533E  + 33 

4.311 

E + 3 3 

1 

6 3 9 E+  3 3 

3 . 53  6E  + 3 2 

2.53-E+-- 

2 . 753E+33 

4.311 

E + 3 3 

1 

61 4E+33 

3.33+E+22 

2. 67  aE+ - . 

3.  333E+33 

4.311 

E + 33 

1 

369E+33 

4.243E+J2 

-. 797E+- 5 

3.  2 52  E+03 

4.311 

E + 33 

i . 

1 23E+33 

4 . j 9 5E  + 52 

2 . 9 5 E + 5 _ 

3. 533E+33 

4.311 

E+33 

3. 

7 3 2 E + 3 3 

-i . • 3 E + 3 _ 

2 . 9 3 2 E+ 3 2 

3 . 7 52  E + 2 <3 

4.311 

E+3  3 

6. 

3 33 E+3 2 

i . ■ 3 E + 

3 . 5 4 4 4 + 5 5 

4. 333E+33 

4.  31  1 

E + 33 

3 • 

3 7 9 E + 32 

. 6 5 7 E + - 2 

4 . - ■ 5 E+ 3 2 

-4.  253  E+33 

4.  3 1 1 

E+33 

1 . 

427E+32 

6.3  1 -E  + 32 

j.1J;E+5- 

4.  53  3 E+3  3 

4.311 

E + 33 

- 1 . 

32  5E+3  2 

6 . 3 64  E + 32 

3 . 1 1 3 E + 3 _ 

4. 7 S3 E + 3 3 

4.311 

E + 33 

- 3 • 

47  6 E+  32 

fc. 71 7E+22 

3.292E+t 2 

5.333E+33 

4.311 

E + 3 3 

- *5. 

9 2 3 E + 3 2 

7.371  E + 3 2 

3 . - 53E+5 . 

5. 2 53E  + 33 

4.311 

E+33 

- 3 . 

3 33  E + 22 

7 . 4 2 5 E + 3 2 

2 . 9 9 j E + . ' _ 

5. 53  3 E + 3 3 

4.311 

E + 33 

- 1 . 

3 3 3 E +3  3 

7. 773E+ 

2 . 9 1 + E + 3 2 

5.  7 53  E + 3 3 

4.311 

E+3  3 

- 1 . 

323E+33 

3 . 1 3 E + : 

2.31 5~+-2 

6.  3 3 3 E «-3  3 

4.311 

E + 33 

- 1 . 

673E  + 3 3 

3.453E+22 

2.  6 9 9 E + 5 2 

6 . 2 53  E+3  3 

4.311 

E+3  3 

- 1 . 

3 1 9 E + 3 3 

_5  • ‘*3  JE+i'  - 

2 . o 6 5 E + 2 - 

6 . 53  3 E +3  3 

4.311 

E + 3 3 

3 64 E +3  3 

9.1  9 2 E + 

2 . 45  1 E + - - 

6.  7 53  £ -t-3  3 

4.711 

E + 1 J 

_ ^ 

33  9E+3  3 

• . 4 E + 

2 . 2 2 2 Z + 

7.  53  ' • * • ' 

4.  3 1 1 

E+3  3 

- 3 • 

5 5 4 E +33 

9 , 39  9E+3  2 

2 . 5 2 2 E + . 

7. 253E+32 

4.  3 1 1 

E + 3 3 

) 

79  9E  + 3 3 

1 . 3 25E+  < 

1 . ; 5 3 E +~  - 

7.  ! ' ’ E *3 3 

4.311 

E + 3 3 

- 3. 

3 4 4 E + 3 3 

1 . 3 fc  1 E + 3 3 

1 . 5 6 4 E + - - 

7.  7 53  E + 3 3 

4.311 

E + 3 3 

- 3. 

293  E + 3 3 

1 . 3 9 6 E + 5 5 

1 . 3 5 4 E+  5 2 

J.333E+33 

4.311 

E + 3 3 

- 3 • 

6 3 3 E + 3 3 

1 . 1 3 1 E + 3 3 

1 . 5_4E+~2 

••  . 2 53 £ + 33 

4.  3 1 1 

E + 33 

- 3 . 

733  E+3  3 

1 . 1 67E+ - 3 

7 . 2472  + . 1 

3.  533 E+33 

4.311 

E+3  3 

- 4. 

32  jE  + 3 3 

1 . 2 - . E + 2 J 

4 . 5 4 5 E + 5 1 

3 . 753E+33 

4.311 

E + 33 

- 4. 

2 73  E + 3 3 

1 . 2 3 7E  + _ j 

6.47.2E+53 
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Trajectory  parameters  at  3 seconds 


Case  1 

Case  2 

K/2 

9.56  x 10~b 

0 

x(  cm) 

1201 1 . 3 

12932. 0 

y(cm) 

7 79  1 .0 

8525.3 

dx/dt( cm/ sec ) 

Cl 

to 

4 311. 0 

dy /dt ( cm/ sec ) 

990.8 

1369.0 

In  all  cases  the  PROSE  program  is  given  an  initial 
estimate  of  the  parameter  sought  and  it  then  iterates 
through  the  equations  of  motion  until  the  desired  testing 
criterion  is  satisfied;  in  this  case  to  minimize  SUMSQR 
in  Equation  23.  As  a first  test,  the  analysis  routine 
was  asked  to  calculate  the  drag  force  in  Case  1 when  all 
other  parameters  were  defined  exactly  (we  ignore  8 initially). 
Table  G is  a listing  of  the  8 iterations  on  the  drag  force 
(called  F)  which  were  needed  to  minimize  SUMSQR  with  an 
initial  guess  of  9.6  x 10  for  F and  a 10  data  point  set. 
Convergence  is  fairly  rapid  (only  25  seconds  of  total 
running  time  on  a remote  teletype  terminal),  the  calculated 
drag  force  is  almost  exact,  and  the  SUMSQR  (or  error)  is 
also  small . Note  if  one  assumes  in  SUMSQR  that  AX  ~ AY, 
then 

N 

SUMSQR  s V 2AX  * 2NAX  (21) 

i = 1 


and 


^ . SUMSQR 
2N 


(25) 


Ta  b 1 e 6 


Case  1 PROSE  Results  for  Drag  Force  Only 


COM VERGEMC  E COMDITIOM  AFTER  3 I TERA  T 10  MS 
UMKMOWMS  COM  VERGED 
03JECTIVE  CR  I TER  10 M SATISFIED 
ALL  SPECIFIED  CRITERIA  SATISFIED 


LOOP  MUMBER 
UMKMOWMS 

F 

OBJECTIVE 
SUM  SO  R 


CINITIAL1  1 2 

9.  633JJJE-J  5 7.6J33J3E-35  5.  63J3JJE-J5 

5.  33967  3E  + 3 7 4.153349E  + 37  2 . 432  57  3 E + 3 7 


LOOP  MUMBER  CIMITIAL1  3 4 

UMKMOWMS 

F 9.63J33JE-J5  3.b33333E-35  1 . 6 33  33  <3  E - 3 5 

OBJECTIVE 

SUMSQR  5.  339673E+3  7 1 . 3 1 1 9 3 5 E+- 3 7 7 . 5 7 7 3 6E  + 3 5 

LOOP  MUMBER  [IMITIAL1  5 6 

UMKMOWMS 

F 9.  633333E-3  5 3.  5331  69E-36  9.539764E-36 

03JECTIVE 

SUMS3R  5.339673E+37  2 . 53 o 3 1 1 E +3 4 1.334333E+33 


LOOP  MUM3ER 
UMKMOWMS 

F 

03JECTI VE 
SUM  S3 R 


CIMITIAL1  7 3 

9.  63  3333  E- 35  9.  565737E-36  9.5b5753E-3b 

5. 3 3 96  73  E +3  7 1.323732E+33  1.323732E+33 


Here  SUMSQR  is  1320.7,  and  with  N = 10,  this  implies  that 
AX  “ 60  cm  over  a portion  ot  the  trajectory  which  extends 
from  12000  cm  to  26000  cm. 

The  PROSE  routine  was  further  exercised  by  asking 
it  next  to  calculate  initial  velocity  components  and  then 
to  add  initial  position  coordinates.  These  results  are 
shown  in  Tables  7 and  S and  it  can  be  seen  that  although 
the  number  of  variables  to  be  determined  increased  from 
1 to  3 and  5,  the  number  of  iterations  required  only  went 
from  8 to  10  and  12.  Note  also  that  in  the  latter  case, 
all  five  variables  have  been  calculated  with  errors  which 
are  considerably  less  than  one  percent.  Note  also  in 
this  case  that  SLMSQR  is  only  4. 15  or  that  the  average 
error  between  the  two  trajectory  coordinates  is  only  0.2 
cm  over  the  same  12000  to  26000  cm  range!  Similar  results 
and  accuracies  were  obtained  for  Case  2 (zero  drag)  as 
is  shown  in  Table1  9,  and  to  test  the  effects  of  the  initial 
estimates  the  5 par  tune  ter  Case1  1 exercise  was  rerun  with 
crude  initial  estimates  as  shown  in  Table  10.  In  the 
latter  case,  convergence  was  accomplished  in  only  9 iterations 
aga i n with  accuracies  better  than  one  percent. 

Having  confirmed  the  validity  of  this  technique 
with  the  above  tests,  trignometric  conversions  to  include 
the  azimuthal  ejection  angle,  3,  were  included  in  the  PROSE 
coding  and  most  of  the  above  tests  were  redone  for  3=0 
to  assure1  that  no  programming  errors  had  been  introduced. 

A new  "ideal"  data  set  was  then  generated  for  8 = 13° 
using  the  two  cases  found  in  Table  5,  and  til  1 subsequent 
testing  seemed  to  indicate  that  extension  to  the  sixth 
variable  could  bo  accomplished,  and  that  similar  accuracies 
could  be  ant icipated. 

However,  it  soon  became  apparent  that  there  is  a 
subtli1  trade-off'  between  azimuthal  ejection  angle  and 
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Case  1 PROSK  Results  for  Drag  and  Velocity 


CONVERGENCE  CONDITION  A FT ER  1 4 ITERATIONS 
'UNKNOWN S CONVERGED 
OBJECTIVE  CRITERION  SATISFIED 
ALL  SPECIFIED  CRITERIA  SATISFIED 


LOO3  N'TMRER  [INITIAL!  1 

UNKNOWNS 

17  9.  64  4 4 44  E- 4 5 - 34  72  57E-4  5 -2 

U4  3. 724444F*4 J A. 471 574E  + 2 3 3 

V4  7.332222E+22  9.744355E+42  1 

OBJECTIVE 

SUN  SO  R 6.  222435E-47  1 . 93 1 2 1 I E <-2  i 2 

LOO3  NUMBER  t INITIAL!  3 

UNKNOWNS 

17  9.644444E-45  - 5 . 0 9 4 5 6 2 E - 4 3 - J 

U4  3.  724424‘7<-43  2.  AS  1 27  s :■>  43 

7 . 3 3 0 2 4 0 E + 0 2 l . 4 6 5 9 36  E * 4 3 I 

OBJECT  I VF 

S'M  33"  6 . 2 I 2 A 3 ^E  ► 4 7 1.422A23E«-4  7 3 


LOO3  NUMBER 
UNKNOWN  S 


U4 

V4 

OBJECTIVE 
SUM  5 OR 


C I N I TI  AU 


9 .643444"-!  5 - 7.  71 974E-46  A, 

3.724444E+43  5.  33  141  E«-33  3. 

5 ■ : * ' 1.44  ' • >43 

6.222435E*47  7.432  1 E+4  7. 


LOO3  NUMBER 
UNKNOWN S 

F 

ua 


OBJECTIVE 
SUMSO  ■ 


: I N I T I A L 1 7 


9.  644444E-45  3. 

3.  7 244 42 E *-4  3 3. 

7.  334444 E+42  •. 


34  3564E-46  9. 

732973E*03  3. 

9 5 34  93  Et- 4 2 9. 


6. 22243  5E*4  7 1.355361 £+43  1. 


LOO3  NUMBER 
UNKNOWN 5 

F 

U4 

U4 

OBJECTIVE 

-OMSOR 


C INITIAL! 

9 , 624444E-4 5 
3. 724324 E+ 43 
7.  334444  Et- 42 

6.2  22  43  E*-4  7 


9 


9.  532  565E-4  6 9. 

1.75  13 1 E +4  3 3 . 

9.9  537 

1 .6463  9 ■ ♦ ' I 1 . 


If. 


3 7 51  3 9 E-  2 5 

4 62  2 4 4 E +■  4 3 
246497E<-23 

5 3 7-.  <3 E *■  4 7 


;a1s44E-45 
9 1 4 9 53  E<-4  3 
4 3 4 6 6 6 E * 4 3 

2 9771  7E*4  6 


6 

4 7 422  9 E-4  fc 
ol 1 7 3 3 " +4  3 
973  1 4 3 F «•  4 2 

4 625  14" *4  4 


51  A 53  5 E-  4 6 
751 597E+ 43 
9 5377aE*42 

67  4 b 34  E *■  4 I 


1 4 

534577E-46 
7 2331 E+43 
953753E+42 

6465  9 6E  t-4  1 


Case  1 PROSE  Results  for  Drag,  Velocity  and  Position 
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Table  9 


Case  2 PROSE  Results  for  Drag,  Velocity  and  Position 


CONVERGENCE  CONDITION  AFTER  3 ITERATIONS 
UNKNOWNS  CONVERGED 
OBJECTIVE  CRITERION  SATISFIED 
ALL  SPECIFIED  CRITERIA  SATISFIED 


LOOP  NUM3ER  

, . . . . C INITIAL] 

1 

2 

UNKNO WN S 

F 

9,  6<3«5  <3^  <3  6 

- 2 . 

2 A3 3 2 1 E-3  S 

1 . 1 6 1 1 26E-2  5 

U3 

a.  31 3223E+23 

A. 

1 1 6336E+3 3 

5.  225665  E *23 

V2 

i . 123333E+33 

1 . 

663352  E + 2 3 

1 . 61  7A76E+3  j 

xa 

1 . 293333E+3A 

1 . 

2 32  1 44E+.2  4 

1 . 22  57  3 6E  + 2 A 

ya 

3. 525232E+23 

7. 

92  665  iE  + 2 3 

3 . 1 476  12  E+2 3 

OBJECTIVE 

SUMS3R 

6.  3231 63  E+  2 6 

1 . 

1 263A2E+2  7 

1 . 323273E+26 

LOO?  NUM3ER  

, C INITIAL] 

3 

A 

UNKNOWNS 

F 

9. 622222 E-26 

2.2551 25E-26 

1 . 2 3572  AE-2 6 

U2 

A. 31 32  22  E + 2 3 

A.  A23733E+23 

A.  3573 A2  E+  2 3 

V2 

1 . 1 23323 E+33 

1 . 322A72E+33 

1 . 37362  AE+2  3 

X2 

1 . 2 93  33  3 E+  3 A 

1 • 2 65  i<3  4E ♦ .3  4 

1 . 23is37E  + 2A 

Y2 

3.  5 2 52 2 2 E+2 3 

:s  . 533967E+33 

3.  525367E+03 
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Case  1 


PROS l-l  Results  with  Crude  Initial  Estimates 
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drag  force  which  makes  it  difficult  to  resolve  the  "correct" 
answer.  This  trade-off  does  not  appear  to  be  continuous, 
but  has  an  apparent  double  value  as  is  typified  by  the 
results  of  Table  11.  The  TROSE  coding  in  this  example 
could  not  optimally  minimize  SL'MSQR,  but  appeared  to  oscil- 
late around  the  correct  answers  calculating  drag  forces 
of  1 x 10-5  and  7.G  x 10- 5 and  azimuthal  ejection  angles 
of  1-1  and  18.6°.  This  situation  is  graphically  demonstrated 
in  Figure  1 1)  where  SU.MSQR  is  plotted  as  a function  of 
iteration  number  and  the  calculated  drag  and  azimuth  are 
also  indicated.  The  coding  is  obviously  not  converging 
on  a single  answer,  and  it  does  not  appear  that  the  situation 
will  improve  even  if  the  iteration  process  is  allowed  to 
proceed  further.  (PROSE  will  automatically  terminate  the 
calcu 1 at ional  process  if  an  answer  has  not  been  obtained 
alter  20  iterations.  It  can  be  overridden . ) Note  also 
that  there  is  a definite  trade-off  in  drag  force  and 
angle  as  even  the  relative1  minimums  for  SL'MSQR  are  hall 
11 0 and  18.6°.  There  are  indications  that  by  appro- 
priately bounding  the  s i ze  of  the  steps  which  the  code 
may  take  in  selected  variables,  that  this  behavior  may 
be  minimized  and  perhaps  even  eliminated,  but  because  of 
time  constraints  we  have  not  been  able  to  verify  this. 

It  might  appear  at  this  point  that  the  goal  of 
obtaining  the  desired  ejection  parameters  from  the  tra- 
jectory data  is  unobtainable.  However,  consider  the 
last  two  iterations  of  Table  11  and  compare  those 
results  with  the  initial  conditions  from  Table  5.  (See 
Table  12. ) Note  here  that  in  iteration  19  the  calculated 
parameters  ol  position  and  velocity  differ  from  the  known 
values  by  a maximum  of  d percent,  while  on  iteration  20 
the  di I ference  is  less  than  one  percent.  What  then  does 
this  imply  with  respect  to  the  dynamic  parameters  at 
the  o r i g i n ? 
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Case  1 PROSE  Result  with  |? 


15°, Showing  Observed  Drag/ Angle  Ambiguity 
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Table  12 


Comparison  of  Known  and  Calculated 
Dynamic  Parameters  at  3 Seconds 


Case  1 (Input) 

I te rat  ion 

19 

Iteration  20 

F 

9.57  x 10  1 ' 

7.00  x 10 

-0 

1.00  x 10“° 

dx/dt ( cm/ sec  ) 

3,753 

3,879 

3,717 

dy /dt ( cm/sec ) 

991 

990 

989 

x ( cm ) 

12,011 

12 , 350 

1 1 ,906 

y ( cm ) 

7,79-1 

7,877 

7,777 

iJ  ( degrees ) 

15 

18.7 

14.3 

The  calculated  parameters  of  Table  12  were  used  as 
initial  conditions  in  the  differential  equations  of  motion 
and  the  respective*  trajectories  were  traced  back  toward  the 
origin.  These  results  are  again  compared  below: 


Table  13 


Comparison  of  Known  and  Calculated 
Dynamic  Parameters  at  0.25  Seconds 


Input 

It.-  ra  t i on  1 9 

I t ora t ion 

20 

F 

9.57  x 10"° 

7.60  x 1 0-( ’ 

1 . 00  x 1 0 

-6 

dx/ dt ( cm/ sec ) 

4249 

4288 

•1229 

dy/dt ( dm/sec  ) 

4006 

3946 

4018 

x ( cm) 

1070 

1 172 

1047 

We  see  here  that  propagating  the  three  percent  errors 
of  Table  12  thru  the  equations  of  motion  result  in  errors 
as  high  as  thirteen  percent,  but  more  generally  of  only 
several  percent.  The  larger  positional  errors  of  iteration 
19  would  place  the  ejecta  at  38  feet  in  x.y  as  opposed  to 
the  correct  3d  to  35  feet,  but  then  ag.in  this  is  probably 
better  accuracy  than  one  could  hope  to  obtain  from  any 
cratering  calculation  (note  that  velocity  errors,  for 
either  iteration,  are  definitely  less  than  two  percent). 

Although  we  have  not  definitively  demonstrated  the 
feasibility  of  obtaining  ejecta  dynamic  parameters  with 
sufficient  accuracy  to  correlate  with  cratering  calculations, 
the  above  demonstration  would  seem  to  imply  that  the 
suggested  technique  has  sufficient  merit  to  warrant  further 
examination.  The  simplicity  of  concept  and  the  lack  of 
preconceived  bounds  or  constraints  on  allowable  solutions 
should  also  prove  appealing.  In  essence,  this  approach 
is  straightforward  curve  fitting  (calculated  versus 
measured  trajectories)  which  utilizes  a 1 1 of  the  acquired 
data  to  generate  a solution  to  a non-linear  optimal 
estimation  problem.  Note  also  that  in  the  generation 
of  SUMSQR  we  automatically  obtain  a measure  of  the  goodness 
of  fit  between  measured  and  calculated  trajectories. 
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CONCLUSIONS  AND  RLCOMMLN DAT IONS 


V . 


The  work  reviewed  in  this  report  has  considered 
in-flight  sizing  for  ejecta  from  the  MIDDLE  GUST  high 
explosive  test  series,  and  has  also  examined  i ti  some 
detail  the  feasibility  of  obtaining  initial  dynamic 
ejecta  characteristics  for  correlation  with  cratering 
calculations.  Both  efforts  have  made  significant  con- 
tributions to  our  understanding  of  the  cratering  process 
and  of  ejecta  dynamics.  Each  is  considered  more  fully 
be  low . 


I n- F 1 i gh  t Sizing 

The  ejecta  sizing  measurements  presented  in  Section 
III  tend  to  support  the  a-'5'0  or  a-1'0  power  law  currently 
employed  in  the  SA 1 ejecta  model.  As  noted  in  this  report 
these  distributions  (and  hence  this  conclusion)  hold  most 
specifically  for  the  farther-out  or  discontinuous  ejecta 
blanket  as  close-in  measurements  could  not  be  made  because 
of  the  dust /combust  ion  products  cloud.  The  support 
which  those  measurements  provide  only  becomes  apparent 
after  detailed  comparison  with  ejecta  fields  predicted 
by  the  SAI  model.  Gravitational  and  aerodynamic  filtering 
of  the  initial  size  distributions  create  the  fields 
which  can  be  used  for  direct  comparisons  such  as  that 
shown  in  figure  1. 

Although  calibration  shot  s) , over  a shale  geology, 
tended  to  produce  more  ejecta  than  the  other  calibration 
shots,  the  general  shape  of  the  distributions  was 
similar.  The  calibration  shots  and  the  20  ton  shot 


were  all  half-buried  and  therefore  the  data  presented  for 
them  are  directly  comparable  as  there  is  no  variation  in 
shot  geometry.  For  these  shots  an  increase  in  yield 
by  a factor  of  40  (0.5  to  20  ton)  increases  e.jecta 
yield  by  a factor  of  10  to  20. 

If  there  is  sufficient  interest  in  the  impact 
characteristics  of  the  discontinuous  ejecta  field,  these 
measurements  should  be  pursued  to  cover  other  geologies 
and  shot  configurations.  At  best,  the  data  acquired  to 
date  does  not  contradict  the  existing  ejecta  size  model, 
but  not  enough  data  has  been  acquired  to  perform  a 
definitive  validation. 

Ejecta  Dynamics 

Although  existing  analysis  techniques  may  provide 
data  on  ejecta  dynamics  which  is  of  sufficient  accuracy 
for  certain  applications,  a better  know  1 edge  of  ejecta 
trajectories  may  provide  data  of  sufficient  accuracy  and 
validity  to  permit  direct  comparison  with  (and  hopefully 
validation  of)  existing  cratering  models.  The  approach 
reviewed  here  has  such  promise  and  does  not  rely  upon 
preconceived  bounds  or  constraints  to  permit  identification 
of  the  "correct"  solution. 

This  problem  has  been  approached  with  a unique 
and  powerful  programming  tool  (PROSE)  which  made  the 
entire  investigation  feasible  within  the  existing  constraints 
Although  this  investigation  has  revealed  an  inherent 
ambiguity  between  the  drag  force  and  the  azimuthal  ejection 
angle,  it  is  an  ambiguity  which  is  not  unique  to  this 
approach.  Other  analysis  techniques  have  the  same  ambiguity, 
but  have  not  been  able  to  recognize  it  because  of  their 
coa rseness . 

ll  is  felt  that  the  approach  suggested  here  has 
suf fit  lent  merit  to  warrant  further  investigation, 


particularly  with  respect  to  developing  tools  which  will 
permit  resolving  this  ambiguity.  There  are  constraints 
^ available,  within  the  PROS!  system,  on  exactly  how  certain 

variables  may  be  approached  and/or  approximated  (not  on 
their  final  value)  and  these  may  offer  a simple  solution 
to  the  problem.  With  the  development  of  a validated 
analysis  technique,  DMA  must  then  decide  whether  or  not 
a correlation  of  ejecta  dynamics  and  cratering  models 
mer  i t s furl  lie r work  . 


REFERENCES 


Seebaugh,  W.R.,  "A  Dynamic  Crater  Ejecta  Model," 
Science  Applications,  Inc.,  SAI-76-654-WA , In 
Press,  October  197(5. 

Layson , W.M.,  et  al . Private  Communication,  General 
Research  Corporation,  August  19(59. 

Gault,  D.E.,  Shoemaker,  E.M.,  and  Moore,  11. T., 

"Spray  Ejected  from  the  Lunar  Surface  by  Meteoroid 
Impact,"  NASA  TND-1767,  April  1963. 

Linnerud,  J.J.,  "In-Flight  Ejecta  Size  Distributions 
for  MIDDLE  GUST  III,"  Science'  Applications,  Inc., 
SAI -74-502-BN , DNA3574F , March  1975. 

Seebaugh,  W.R.,  "Studies  of  the  Nuclear  Crater 
Ejecta  Environment,"  Science  Applications,  Inc.. 

S A I - 7 5- 50  7- W A , DNA3601F,  June  1975. 

Wisotski,  J.,  "MIDDLE  GUST  Dynamic  Ejecta  Measure- 
ments," Denver  Research  Institute,  In  Press , 

July  1974. 

Linnerud,  11. J.,  "Ejecta  Sizing  for  MIDDLE  GIST  III 
--  A Feasibility  Demonstration,"  Science  Appli- 
cations, Inc..  SAI-73-503-BN , July  1973. 

Linnerud,  '1.  J..  "MIDDLE  GL'S’l  1 Eject  a Sizing  and 
an  Examination  « > t Ejecta  Ballistics,"  Science 
Applications  Inc.  , SA  1-7(5-504-150 . DNA  Final  Report  , 
In  Press,  February  197(5. 

Wisotski,  J..  "MIDDLE  NORTH  Series.  MIXED  COMPANY 
Events,"  Denver  Research  Institute,  DRI2659, 

POR6625 , August  1975. 


APPENDIX 


TABULATED  SIZING  DATA  FOR  MG  CALIBRATION  EVENTS 

TABLE  1 

INTEGRAL  EJECTA  COUNT  FOR  CALIBRATION  SHOT  3a 


Area 

o 

PP~ 


2 

cm 


2 sec 
West 


4 sec 
West 


10 

1 . 3 

3722 

6131 

20 

2.6 

3271 

5030 

40 

5 . 3 

2787 

3781 

60 

7.9 

2477 

3024 

80 

10.6 

2237 

2473 

100 

13.2 

2022 

2088 

200 

26.4 

1378 

1149 

300 

39.6 

1011 

774 

400 

52 . 8 

780 

549 

800 

105 . 6 

370 

216 

1000 

132.0 

301 

151 

2000 

263.9 

104 

36 

59 


A 


m 


TABLE  2 


L 


INTEGRAL  EJECTA  COUNT  FOR  CALIBRATION  SHOT  (5 


Area 


PP 


cm 


sec  £ sec 


West 


Eas  t 


4 sec 
West 


6 sec 
West 


table 


b 


I NTEGRAL 

EJECT A COUNT 

FOR  CALIBRATION 

SHOT 

8 

Area 

2 sec 
West 

2 sec  4 

East 

sec 
Wes  t 

6 sec 
West 

2 

PP 

2 

cm 

10 

3.7 

1751 

2002 

2528 

2776 

20 

7 . 5 

1503 

1668 

2146 

2 196 

■10 

15.0 

1319 

1343 

1764 

1498 

60 

22.5 

1202 

1137 

1517 

1172 

HO 

30.0 

1101 

995 

1 339 

952 

100 

37.6 

1021 

869 

1 182 

817 

200 

75 . 1 

711 

52  3 

692 

4 34 

300 

112 .6 

550 

348 

482 

269 

•100 

150.2 

160 

253 

333 

197 

800 

300.1 

2 17 

9 1 

109 

58 

1000 

375.5 

190 

69 

56 

40 

2000 

75  1 . 0 

88 

19 

8 

2 

6 1 


TABLE  1 


INTEGRAL 

EJECTA  COUNT 

FOR  CALIBRATION  SHOT 

9 

Area 

2 sec 
West 

2 sec 
East 

4 sec 
West 

6 sec 
West 

2 

PP 

2 

cm 

10 

3.8 

6612 

6506 

5209 

3025 

20 

7 . 6 

4973 

5 196 

3996 

1915 

40 

15.1 

3659 

3969 

294  1 

1218 

60 

22.7 

3012 

3215 

2408 

885 

80 

30 . 3 

2534 

2682 

2044 

683 

100 

37.9 

2218 

2299 

1764 

533 

200 

75.7 

1310 

1283 

983 

220 

300 

113.6 

894 

802 

570 

124 

400 

151. 4 

64  0 

569 

575 

79 

800 

302.9 

274 

194 

89 

22 

1000 

378.6 

1 75 

145 

49 

12 

2000 

757.2 

38 

32 

1 

3 
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